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ABSTRACT 

Two-dimensional (2D) hydrodynamical simulations of the deleptonization of a newly 
formed neutron star (NS) were performed. Driven by negative lepton fraction and en- 
tropy gradients, convection starts near the neutrinosphere about 20-30 ms after core 
bounce, but moves deeper into the protoneutron star (PNS), and after about one sec- 
ond the whole PNS is convective. The deleptonization of the star proceeds much faster 
than in the corresponding spherically symmetrical model because the lepton flux and 
the neutrino (u) luminosities increase by up to a factor of two. The convection below 
the neutrinosphere raises the neutrinospheric temperatures and mean energies of the 
emitted z/'s by 10-20%. This can have important implications for the supernova (SN) 
explosion mechanism and changes the detectable v signal from the Kelvin-Helmholtz 



cooling of the PNS. In particular, the enhanced v e flux relative to the v e flux dur- 
ing the early post-bounce evolution might solve the overproduction problem of certain 
elements in the neutrino-heated ejecta in models of type-II SN explosions. 

Subject headings: supernovae: general - stars: neutron - elementary particles: neutri- 
nos - turbulence - convection - hydrodynamics 

1 Introduction 

Convection in the newly formed NS might play an important role to explain the ex- 
plosion of a massive star in a type-II SN. Epstein (1979) pointed out that not only en- 
tropy, S, inversions but also zones in the post-collapse core where the lepton fraction, 
Yi, decreases with increasing radius tend to be unstable against Ledoux convection. 
Negative S and/or Y t gradients in the neutrinospheric region and in the layers between 
the nascent NS and the weakening prompt shock front were realized in a variety of 
post-bounce SN models by Burrows & Lattimer (1988), and after shock stagnation in 
computations by Hillebrandt (1987) and more recently by Bruenn (1993), Bruenn & 
Mezzacappa (1994), and Bruenn et al. (1995). Despite different equations of states 
(EOS), v opacities, and v transport methods, the development of negative Y\ and S 
gradients is common in these simulations and can also be found in PNS cooling models 
of Burrows & Lattimer (1987), Keil & Janka (1995), and Sumiyoshi et al. (1995). 

Convection above the neutrinosphere but below the neutrino-heated region can 
hardly be a direct help for the explosion (Bethe et al. 1987, Bruenn et al. 1995), whereas 
convectively enhanced lepton number and energy transport inside the neutrinosphere 
raise the v luminosities and can definitely support neutrino-energized SN explosions 
(Bethe et al. 1987). In this context, Burrows (1987) and Burrows & Lattimer (1988) 
have discussed entropy-driven convection in the PNS on the basis of ID, general rel- 
ativists (GR) simulations of the first second of the evolution of a hot, 1.4 M Q PNS. 
Their calculations were done with a Henyey-like code using a mixing-length scheme for 
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convective energy and lepton transport. Recent 2D models (Herant et al. 1994, Bur- 
rows et al. 1995, Janka & Miiller 1996 and references therein) confirmed the possibility 
that convective processes can occur in the surface region of the PNS immediately after 
shock stagnation ("prompt convection") for a period of at least several 10 ms. These 
models, however, have been evolved only over rather short times or with insufficient 
numerical resolution in the PNS or with a spherically symmetrical description of the 
core of the PNS that was in some cases even replaced by an inner boundary condition. 

Mayle & Wilson (1988) and Wilson & Mayle (1988, 1993) demonstrated that con- 
vection in the nascent NS can be a crucial ingredient that leads to successful delayed 
explosions. With the high-density EOS and treatment of the v transport used by the 
Livermore group, however, negative gradients of Y\ tend to be stabilized by positive S 
gradients (see, e.g., Wilson & Mayle 1989). Therefore they claim doubly diffusive neu- 
tron finger convection to be more important than Ledoux convection. Doubts about 
the presence of doubly diffusive instabilities, on the other hand, were recently raised 
by Bruenn & Dineva (1996). Bruenn & Mezzacappa (1994) and Bruenn et al. (1995) 
also come to a negative conclusion about the relevance of prompt convection in the 
neutrinospheric region. Although their post-bounce models show unstable S and Y\ 
stratifications, the mixing-length approach in their ID simulations predicts convective 
activity inside and around the neutrinosphere to be present only for 10-30 ms after 
bounce and to have no significant impact on the v fluxes and spectra when an elaborate 
multi-group flux-limited diffusion method is used for the v transfer. Such conclusions 
seem to be supported by recent 2D simulations of Mezzacappa et al. (1996). These 2D 
models, however, still suffer from the use of an inner boundary condition at a fixed ra- 
dius of 20-30 km, a simplified treatment of neutrino-matter interactions, and imposed 
neutrino fluxes and spectra from spherically symmetrical models. 

From these differing and partly contradictory results it is evident that the question 
whether, where, when, and how long convection occurs below the neutrinosphere seems 
to be a matter of the EOS, of the core structure of the progenitor star, of the shock 
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properties and propagation, and of the v opacities and the v transport description. In 
this Letter we compare ID simulations with the first self-consistent 2D models that 
follow the evolution of the newly formed NS for more than a second, taking into account 
the GR gravitational potential and making use of a flux-limited equilibrium diffusion 
scheme that describes the transport of v e , v e , and v x (sum of v^, z/ M , u T , and 9 T ) and is 
very good at high optical depths but only approximative near the PNS surface (Keil 
& Janka 1995). Our simulations demonstrate that Ledoux convection may continue in 
the PNS for a long time and can involve the whole star after about one second. 

2 Numerical implementation 

The simulations were performed with the explicit Eulerian hydrodynamics code Prometheus 
(Fryxell, Miiller, & Arnett 1989) that employs a Riemann-solver and is based on the 
Piecewise Parabolic Method (PPM) of Colella & Woodward (1984). A moving grid with 
100 nonequidistant radial zones (initial outer radius ~ 60 km, final radius ~ 20 km) 
and with up to 60 angular zones was used, corresponding to a radial resolution of a 
few 100 m (^ 1 km near the center) and a maximum angular resolution of 1.5°. In the 
angular direction, periodic boundary conditions were imposed at ±45° above and below 
the equatorial plane. The stellar surface was treated as an open boundary where the 
velocity was calculated from the velocity in the outermost grid zone, the density profile 
was extrapolated according to a time-variable power law, and the corresponding pres- 
sure was determined from the condition of hydrostatic equilibrium. The Prometheus 
code was extended for the use of different time steps and angular resolutions in dif- 
ferent regions of the star. Due to the extremely restrictive Courant-Friedrichs-Lewy 
(CFL) condition for the hydrodynamics, the implicit v transport was computed typ- 
ically with 10 times larger time steps than the smallest hydrodynamics time step on 
the grid (~ 10~ 7 s) (Keil 1996). 

Our simulations were started with the ~ 1.1 M (baryonic mass) central, dense 
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part (p <: 10 n g/cm 3 ) of the collapsed core of a 15 M progenitor star (Woosley et 
al. 1988) that was computed to a time of about 25 ms after core bounce (i.e., a few 
ms after the stagnation of the prompt shock) by Bruenn (1993). Accretion was not 
considered but additional matter could be advected onto the grid through the open 
outer boundary. In the 2D run, Newtonian asphericity corrections were added to the 
spherically symmetrical GR gravitational potential: $ 2 d = ®yd + (^zd ~~ ^id) • This 
should be a sufficiently good approximation because convective motions produce only 
local and minor deviations of the mass distribution from spherical symmetry. Using the 
GR potential ensured that transients due to the mapping of Bruenn's (1993) relativistic 
ID results to our code were very small. When starting our 2D simulation, the radial 
velocity (under conservation of the local specific total energy) was randomly perturbed 
in the whole PNS with an amplitude of 0.1%. The thermodynamics of the NS medium 
was described by the EOS of Lattimer & Swesty (1991) which yields a physically 
reasonable description of nuclear matter below about twice nuclear density and is thus 
suitable to describe the interior of the considered low-mass NS (M ns ^ 1.2 M ). 

The v transport was carried out in radial direction for every angular zone of the 
finest angular grid. Angular transport of neutrinos was neglected. This underestimates 
the ability of moving buoyant fluid elements to exchange lepton number and energy 
with their surroundings and is only correct if radial radiative and convective transport 
are faster. Moreover, v shear viscosity was disregarded. For z/-n, p scattering and 3 
flavors of nondegenerate v and v in local thermal equilibrium with the matter one 
estimates at low densities (p ^ 10 14 g/cm 3 ) a dynamic shear viscosity of rf^ A ~ 1.2 • 
10 22 T 1 2 /p 1 4gcm _1 s _1 (van den Horn & van Weert 1981), and in degenerate nuclear 
matter if v ~ 6.7- 10 22 f(Y p )T l0 /p{{ 3 gem 1 s 1 (Thompson & Duncan 1993) when T 10 = 
77(10 MeV) and p u = p/(10 14 g/cm 3 ). The expression f(Y p ) « 0.63... 1 is a function 
of the proton fraction Y p . With convective velocities v c ~ 10 8 cm/s and length scales 
of convective mixing l c ~ 10 5 cm, one obtains Reynolds numbers H„ = v c l c p/r] u ~ 
10 4 ... 10 5 for typical temperatures and densities in the PNS, and corresponding viscous 
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damping timescales of convective motions of the order of t, 



1/ 



'V 



5-50 s. 



Neutrino viscosity is not extremely important on timescales of ~ 1 s and should not be 
able to suppress turbulent convection (Thompson Sz Duncan 1993, Burrows & Lattimer 



The effective numerical viscosity of the PPM scheme is estimated from Eq. (3.6) 
in Porter & Woodward (1994) to be rj Q / p ~ v c l c 10 y / (27r) 2 , which depends on the 
grid resolution via the parameter y defined as the ordinate in Fig. 1 of Porter & 
Woodward (1994). For structures of size Z c /(10 5 cm) = Z c .5 ^ 1 that are typically re- 
solved by about 10 zones, and for flows with advective Courant numbers C a 0.08, 
Fig. 1 of Porter & Woodward (1994) leads to a dynamic shear viscosity of i] n ^ 
5 • 10 23 pi4/ Ci 5W Ci 8 g cm^'s -1 when f Cj § = f c /(10 8 cm/s). The corresponding Reynolds 
numbers are TZ n ~ (2vr) 2 10^ £ 2000, and Eq. (3.5) of Porter & Woodward (1994) 
yields for the viscous damping timescale t n ~ ~ • 10~ y l c /v c <: 25 / C)5 /f c>8 ms. This 
means that small structures represented by only a few grid cells will be affected by the 
numerical viscosity, but damping timescales for fluid motions on scales l c ^ 10 5 cm are 
still a factor | • 10~ y ^ 25 longer than the overturn timescales t Q ~ l c /v c ~ few ms. 



Convection can be driven by a radial gradient of the entropy per nucleon 5* and/or 
by a gradient of the lepton number per baryon Y\ (Epstein 1979) where Y\ includes 
contributions from e~ and e + and from v e and v e if the latter are in equilibrium with 
the matter. Convective instability in the Ledoux approximation sets in when 



Initially, this criterion is fulfilled between ~ 0.7 M & and ~ 1.1 M Q (black area in Fig. |l|; 
see also Bruenn et al. 1995) and convective activity develops within ~ 10 ms after the 
start of the 2D simulation. About 30 ms later the outer layers become convectively 
stable which is in agreement with Bruenn & Mezzacappa (1994). In our 2D simulation, 



1988). 



3 Results 




(i) 
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however, the convectively unstable region retreats to mass shells ^ 0.9 M & and its inner 
edge moves deeper into the neutrino-opaque interior of the star, following a steeply 
negative lepton gradient that is advanced towards the stellar center by the convectively 
enhanced deleptonization of the outer layers (Figs, p] and |2]). Note that the black area 
in Fig. [I] and the thick solid lines in Fig. |2] mark not only those regions in the star which 
are convectively unstable but also those which are only marginally stable according to 
the Ledoux criterion of Eq. (|l|) for angle-averaged S and Yi, i.e., regions where Cl(V) > 
a ■ max r (|CL(r)|) with a = 0.05 holds. For a ^ 0.1 the accepted region varies only 
little with a and is always embedded by the grey-shaded area where \vg\ > 10 7 cm/s. 
Yet, only sporadically and randomly appearing patches in the convective layer fulfill 
Eq. ([!]) rigorously. Figure shows that the black region in Fig. [I] coincides with the 
layers where convective mixing flattens the S and Y\ gradients. 

The convective pattern is extremely non-stationary and has most activity on large 
scales with radial coherence lengths of several km up to ~ 10 km and convective "cells" 
of 20°-30° angular diameter, at some times even 45° (Fig. |3|). Significant over- and 
undershooting takes place (grey regions in Fig. |l|) and the convective mass motions 
create pressure waves and perturbations in the convectively stable NS interior and in 
the surface layers. The maximum convective velocities are usually ~ 4 • 10 8 cm/s, but 
peak values of ~ 10 9 cm/s can be reached. These velocities are typically 5-10% of 
the average sound speed in the star. The kinetic energy of the convection is several 
10 49 erg at t ^ Is and climbs to ~ 2 • 10 50 erg when the PNS is fully convective. 
Relative deviations of Y/ from the angular mean can be several 10% (even 100%) in 
rising or sinking buoyant elements, and for S can reach 5% or more. Rising flows always 
have larger Yi and S than their surroundings. Corresponding temperature and density 
fluctuations are only ~ 1-3%. Due to these properties and the problems in applying 
the Ledoux criterion with angle-averaged S and Y\ straightforwardly, we suspect that 
it is hardly possible to describe the convective activity with a mixing-length treatment 
in a ID simulation. 
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Our 2D simulation shows that convection in the PNS can encompass the whole 
star within ~ 1 s and can continue for at least as long as the deleptonization takes 
place, possibly even longer. A deleptonization "wave" associated with the convectively 
enhanced transport moves towards the center of the PNS. This reduces the timescale for 
the electron fraction Y e to approach its minimum central value of about 0.1 from ~ 10 s 
in the ID case, where the lepton loss proceeds much more gradually and coherently, to 
only ~ 1.2 s in 2D. With convection the entropy and temperature near the center rise 
correspondingly faster despite a similar contraction of the star in ID and 2D (Fig. |]). 
Convection increases the total lepton number flux and the v luminosities by up to a 
factor of 2 (Fig. [5]) and therefore the emitted lepton number iVj and energy E u rise 
much more rapidly (Fig. ^). The convective energy (enthalpy plus kinetic energy) flux 
dominates the diffusive v energy flux in the convective mantle after t ^ 250 ms and 
becomes more than twice as large later. Since convection takes place somewhat below 
the surface, z/'s take over the energy transport exterior to ~ 0.9 M & . Thus the surface 
v flux shows relative anisotropies of only 3-4%, in peaks up to ~ 10%, on angular 
scales of 10°-40°. Averaged over all directions, the neutrinospheric temperatures and 
mean energies (e Ui ) of the emitted v e and v e are higher by 10-20% (Fig. |). 

4 Consequences and conclusions 

The increase of the z/ e -luminosity relative to the z/ e -luminosity during t ^ 0.4 s (Fig. |||) 
will raise in the neutrino-heated SN ejecta. If weak equilibrium is established, 
a particles are absent, and captures can be ignored, captures of v e on n and v e 
on p determine F e ej « 1/ [1 + (L Pe (e Pe ) )/(£„. <£„„))] (Qian & Woosley 1996). The lu- 
minosity ratio enters crucially, since (ep e )/(e^ e ) ~ T Pe /T Ue does not change much due 
to convection (compare Fig. fty. The latter fact can be understood by an analytical 
neutrino Eddington atmosphere model (Schinder & Shapiro 1982) which yields for 
the temperatures T u . of the z/j energyspheres as functions of the effective temperature 



T eff : T v . = ^J2 T cS » 4.6 ^X 52 /r 6 2 MeV, where £ 52 = (L„ e + L P J/(10 52 erg/s), 
r 6 = r/(10 6 cm). The expression £ u . depends on the p and n abundance fractions, Y p and 
Y n « 1 - Yp, in the NS atmosphere. For i/ c one gets ^ e «l + l/(Y ny Jl + 0.182^/^) 
and for P c one finds & e « 1 + 1/(^^1 + 0.210^/^). The ratio T Pe /T, e = ^W^ e 
varies only weakly with the atmospheric composition. Moreover, one can derive that 

increases if / n > jf 4 for f c = £ 2 d/£id and /„ = Md/Md with A/" = L u J(e Ue ) - 
L Pe /(e Pe ). This is fulfilled at times i ^ 0.4 s. The expected increase of Y e ej might help 
to avoid the overproduction of N = 50 nuclei in current SN models (see Hoffman et 
al. 1996, McLaughlin et al. 1996). At times t ^ 1 s the accelerated neutronization of 
the PNS will lead to a more rapid increase of (e Pe ) relative to (e„ e ) than in ID models. 
This will favor a faster drop of Y^ and thus the n-rich conditions required for a possible 
r-processing in the neutrino-driven wind. 

Future simulations of Ledoux convection in the PNS that include the progenitor 
star outside the nascent NS will have to reveal the effects on the SN explosion mech- 
anism. PNS models with (baryonic) masses M ns 1.5... 1.65 M must be considered 
to investigate implications for the v signal detected from SN 1987A. A more elaborate 
description of the v transport and the use of different EOSs are also required. Convec- 
tion in the PNS influences the structure of NS magnetic fields (Thompson & Duncan 
1993), produces gravitational wave emission, and can cause NS recoils by anisotropic 
v emission (Thompson & Duncan 1993). For non-stationary convection with typical 
coherence lengths l c ~ 10 5 cm and overturn timescales t Q ^ R n s/v c & 10 ms, one esti- 
mates a stochastic anisotropy of a ~ (/ c /i? ns )^/t /t ns ^ 10~ 2 (t QS : NS cooling timescale) 
which leads to kick velocities i> ns ~ aE u /(M ns c) of a few lOOkm/s, dependent on the 
energy E v ^ |GM 2 s /i? ns emitted anisotropically in neutrinos. 
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Figure captions: 



FIG. 1. Convective (baryon) mass region inside the PNS vs. time for the 2D simu- 
lation. Black indicates regions which are Ledoux unstable or only marginally stable, 
grey denotes over- and undershooting regions where \vg\ > 10 7 cm/s. 

FIG. 2. Angle- averaged S and Y\ profiles in the PNS. Thick solid lines indicate regions 
that are unstable or only marginally stable against Ledoux convection, crosses mark 
boundaries of over- and undershooting regions where \v e \ > 10 7 cm/s. 

FIG. 3. Panels a and b show the absolute values of the velocity for the 2D simulation 
at times t = 0.525 s and t = 1.047 s, respectively, color-coded in units of 10 8 cm/s. 
The computation was performed in an angular wedge of 90° between +45° and —45° 
around the equatorial plane. The PNS has contracted to a radius of about 21 km at the 
given times. Panels c and d display the relative deviations of the electron fraction Y e 
from the angular means (Y e ) at each radius for the same two instants. The maximum 
deviations are of the order of 30%. Lepton-rich matter rises while deleptonized material 
sinks in. Comparison of both times shows that the inner edge of the convective layer 
moves inward from about 8.5 km at t = 0.525 s to less than 2 km at t = 1.047 s. 

FIG. 4. Radius of the M = 1 M Q mass shell and total lepton number N\ and energy 
E v radiated away by z/s vs. time for the 2D (solid) and ID (dotted) simulations. 

FIG. 5. u e and u e luminosities and mean energies vs. time for the 2D simulation (solid) 
compared with the ID run (dotted). 
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